x=2
y=6
2*x-y
[1] -2
T=c(160, 180, 165, 190, 170, 140, 155, 167, 170, 169)
P=c(80, 75, 70, 85, 83, 87, 70, 89, 75, 85)
S=c("H","F", "F", "H", "H", "F", "F", "F", "H", "H")
C=c("bleue","noir","gris","bleue", "bleue", "gris", "bleue", "noir", "noir", "gris")
df_mat=cbind(T,P)
df_mat
T P
[1,] 160 80
[2,] 180 75
[3,] 165 70
[4,] 190 85
[5,] 170 83
[6,] 140 87
[7,] 155 70
[8,] 167 89
[9,] 170 75
[10,] 169 85
#NB: Les element d'une matrice sont de même type
cbind(df_mat, S, C)
T P S C
[1,] "160" "80" "H" "bleue"
[2,] "180" "75" "F" "noir"
[3,] "165" "70" "F" "gris"
[4,] "190" "85" "H" "bleue"
[5,] "170" "83" "H" "bleue"
[6,] "140" "87" "F" "gris"
[7,] "155" "70" "F" "bleue"
[8,] "167" "89" "F" "noir"
[9,] "170" "75" "H" "noir"
[10,] "169" "85" "H" "gris"
160
180
165
190
170
140
155
167
170
169
80
75
70
85
83
87
70
89
75
85
H
F
F
H
H
F
F
F
H
H
bleue
noir
gris
bleue
bleue
gris
bleue
noir
noir
gris
df=data.frame(T, P, S, C)
head(df)
#Lecture des données dans le fichier DBM1.csv
data=read.csv("Ressources/DBM1.csv", sep = ";")
#Ou utiliser la commande et ne plus convertir en facteur les variables categorielles
data=read.csv("Ressources/DBM1.csv", sep=";", stringsAsFactors = TRUE)
#Conversion en facteur des variables categorielles
data$S=as.factor(data$S)
data$C=as.factor(data$C)
head(data)
mean(data$T)
[1] 174.3111
mean(data$P)
[1] 65.08889
median(data$T)
[1] 177
median(data$F)
[1] 2
quantile(data$T)
0% 25% 50% 75% 100%
157 167 177 181 190
quantile(data$F)
0% 25% 50% 75% 100%
0 1 2 3 7
quantile(data$P, probs = 0.8)
80%
72
min(data$F)
[1] 0
max(data$P)
[1] 98
library(conflicted)
#install.packages("DescTools")
library(DescTools)
Warning: le package 'DescTools' a été compilé avec la version R 4.3.3
Mode(data$S, na.rm = TRUE)
[1] h
attr(,"freq")
[1] 30
Levels: f h
Mode(data$C)
[1] brun
attr(,"freq")
[1] 23
Levels: bleu brun gris noir vert
summary(data)
T P S F C
Min. :157.0 Min. :47.00 f:15 Min. :0.000 bleu:12
1st Qu.:167.0 1st Qu.:59.00 h:30 1st Qu.:1.000 brun:23
Median :177.0 Median :65.00 Median :2.000 gris: 1
Mean :174.3 Mean :65.09 Mean :1.933 noir: 2
3rd Qu.:181.0 3rd Qu.:71.00 3rd Qu.:3.000 vert: 7
Max. :190.0 Max. :98.00 Max. :7.000
#install.packages("psych")
library(psych)
Warning: le package 'psych' a été compilé avec la version R 4.3.3
describe(data[,c("T", "P", "F")])
Desc(data)
------------------------------------------------------------------------------
Describe data (data.frame):
data frame: 45 obs. of 5 variables
45 complete cases (100.0%)
Nr ColName Class NAs Levels
1 T integer .
2 P integer .
3 S factor . (2): 1-f, 2-h
4 F integer .
5 C factor . (5): 1-bleu, 2-brun, 3-gris, 4-noir, 5-vert
------------------------------------------------------------------------------
1 - T (integer)
length n NAs unique 0s mean meanCI'
45 45 0 26 0 174.31 171.66
100.0% 0.0% 0.0% 176.97
.05 .10 .25 median .75 .90 .95
160.20 162.40 167.00 177.00 181.00 184.00 185.00
range sd vcoef mad IQR skew kurt
33.00 8.84 0.05 10.38 14.00 -0.26 -1.20
lowest : 157, 158, 160, 161, 162
highest: 183 (4), 184 (2), 185 (2), 188, 190
heap(?): remarkable frequency (15.6%) for the mode(s) (= 180)
' 95%-CI (classic)
------------------------------------------------------------------------------
2 - P (integer)
length n NAs unique 0s mean meanCI'
45 45 0 21 0 65.09 62.40
100.0% 0.0% 0.0% 67.78
.05 .10 .25 median .75 .90 .95
52.00 54.60 59.00 65.00 71.00 74.20 75.00
range sd vcoef mad IQR skew kurt
51.00 8.95 0.14 8.90 12.00 0.74 2.27
lowest : 47, 49, 52 (2), 53, 57 (3)
highest: 72 (5), 73, 75 (3), 78, 98
heap(?): remarkable frequency (13.3%) for the mode(s) (= 65)
' 95%-CI (classic)
------------------------------------------------------------------------------
3 - S (factor - dichotomous)
length n NAs unique
45 45 0 2
100.0% 0.0%
freq perc lci.95 uci.95'
f 15 33.3% 21.4% 47.9%
h 30 66.7% 52.1% 78.6%
' 95%-CI (Wilson)
------------------------------------------------------------------------------
4 - F (integer)
length n NAs unique 0s mean meanCI'
45 45 0 8 8 1.93 1.45
100.0% 0.0% 17.8% 2.41
.05 .10 .25 median .75 .90 .95
0.00 0.00 1.00 2.00 3.00 4.00 4.80
range sd vcoef mad IQR skew kurt
7.00 1.60 0.83 1.48 2.00 1.05 1.08
value freq perc cumfreq cumperc
1 0 8 17.8% 8 17.8%
2 1 12 26.7% 20 44.4%
3 2 12 26.7% 32 71.1%
4 3 7 15.6% 39 86.7%
5 4 3 6.7% 42 93.3%
6 5 1 2.2% 43 95.6%
7 6 1 2.2% 44 97.8%
8 7 1 2.2% 45 100.0%
' 95%-CI (classic)
------------------------------------------------------------------------------
5 - C (factor)
length n NAs unique levels dupes
45 45 0 5 5 y
100.0% 0.0%
level freq perc cumfreq cumperc
1 brun 23 51.1% 23 51.1%
2 bleu 12 26.7% 35 77.8%
3 vert 7 15.6% 42 93.3%
4 noir 2 4.4% 44 97.8%
5 gris 1 2.2% 45 100.0%
var(data$T)
[1] 78.08283
var(data$P)
[1] 80.17374
var(data$F)
[1] 2.563636
sd(data$T)
[1] 8.836449
sd(data$P)
[1] 8.953979
sd(data$F)
[1] 1.601136
max(data$T)-min(data$T)
[1] 33
max(data$P)-min(data$P)
[1] 51
max(data$F)-min(data$F)
[1] 7
#install.packages("ggplot2")
library(ggplot2)
Warning: le package 'ggplot2' a été compilé avec la version R 4.3.3
ggplot(data, aes(x=F))+geom_histogram(bins = length(unique(data$F)), col="black", fill = "lightgray")+
labs(x = "Nombre de frere", y = "Frequence", title="Distribution suivant le nombre de freres")
ggplot(data, aes(x=F))+geom_boxplot( fill = "lightgray")+
labs(x = "Nombre de freres", title="Boxplot du nombre de freres des individus")
ggplot(data, aes(x=P))+geom_boxplot( fill = "lightgray")+
labs(x = "Poids", title="Boxplot du poids des individus")
ggplot(data, aes(x=T))+geom_boxplot( fill = "lightgray")+
labs(x = "Taille", title="Boxplot de la taille des individus")
ggplot(data, aes(x=T))+geom_histogram(aes(y=..density..),fill="lightgray", col="black")+
geom_density(linewidth=0.55, col="black")+
labs(x = "Taille", y = "Frequence", title="Distribution suivant la taille")
Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data, aes(x=S))+geom_bar(fill="chocolate")+
labs(x = "Sexe", y = "Frequence", title="Distribution suivant le sexe")
tab_C=table(data$C)
#On recupere les lables
labs_C=rownames(tab_C)
labs_C=paste(labs_C, round(100*tab_C/sum(tab_C), 2), sep=" (")
labs_C=paste(labs_C, "", sep="%)")
pie(tab_C,
col=c("#999999", "#E69F00", "#56B4E9","#E9967A","#FFA07A"),
main="Distribution de la couleur des yeux")
legend("topright", labs_C, cex=0.8,fill=rainbow(length(tab_C)))
#install.packages("plotrix")
library(plotrix)
pie3D(tab_C, explode=0.08,
border = "white",
height = 0.25
,main="Distribution de la couleur des yeux")
legend("topright", labs_C, cex=0.5,fill=rainbow(length(tab_C)),
ncol = 1)
#Covariance avec la methode de pearson
print("pearson")
[1] "pearson"
cov(data[,c("T","P","F")], method ="pearson" )
T P F
T 78.082828 45.698990 -2.001515
P 45.698990 80.173737 -4.334848
F -2.001515 -4.334848 2.563636
#Covariance avec la methode de kendall
print("kendall")
[1] "kendall"
cov(data[,c("T","P","F")], method ="kendall" )
T P F
T 1902 764 -214
P 764 1888 -392
F -214 -392 1612
#Covariance avec la methode de spearman
print("spearman")
[1] "spearman"
cov(data[,c("T","P","F")], method ="spearman" )
T P F
T 171.59091 104.5227 -27.46591
P 104.52273 171.5568 -51.75000
F -27.46591 -51.7500 164.36364
#Correlation avec la methode de pearson
print("pearson")
[1] "pearson"
cor(data[,c("T","P","F")], method ="pearson" )
T P F
T 1.0000000 0.5775808 -0.1414663
P 0.5775808 1.0000000 -0.3023637
F -0.1414663 -0.3023637 1.0000000
#Correlation avec la methode de kendall
print("kendall")
[1] "kendall"
cor(data[,c("T","P","F")], method ="kendall" )
T P F
T 1.0000000 0.4031690 -0.1222154
P 0.4031690 1.0000000 -0.2246997
F -0.1222154 -0.2246997 1.0000000
#Correlation avec la methode de spearman
print("spearman")
[1] "spearman"
cor(data[,c("T","P","F")], method ="spearman" )
T P F
T 1.0000000 0.6091996 -0.1635475
P 0.6091996 1.0000000 -0.3081793
F -0.1635475 -0.3081793 1.0000000
library(heatmaply)
Le chargement a nécessité le package : plotly
Warning: le package 'plotly' a été compilé avec la version R 4.3.3
Attachement du package : 'plotly'
L'objet suivant est masqué depuis 'package:ggplot2':
last_plot
L'objet suivant est masqué depuis 'package:stats':
filter
L'objet suivant est masqué depuis 'package:graphics':
layout
Le chargement a nécessité le package : viridis
Warning: le package 'viridis' a été compilé avec la version R 4.3.3
Le chargement a nécessité le package : viridisLite
Registered S3 methods overwritten by 'registry':
method from
print.registry_field proxy
print.registry_entry proxy
======================
Welcome to heatmaply version 1.5.0
Type citation('heatmaply') for how to cite the package.
Type ?heatmaply for the main documentation.
The github page is: https://github.com/talgalili/heatmaply/
Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
You may ask questions at stackoverflow, use the r and heatmaply tags:
https://stackoverflow.com/questions/tagged/heatmaply
======================
#Visualisation de la matrice de correlation
heatmaply(cor(data[,c("T","P","F")], method ="spearman" ),
main = "Matrice de correlation")
#Visualisation de la matrice de covariance
heatmaply(cov(data[,c("T","P","F")], method ="spearman" ),
main = "Matrice de covariance")
#Nuage de points
model=lm(T~P, data)
coefs=coef(model)
ggplot(data, aes(x=P, y=T))+geom_point(col="blue")+
geom_abline(intercept = coefs[1], slope=coefs[2], linewidth=1)+
labs(title =)+
labs(title="Taille en fonction du poids", x="Taille",
y="Poids")
#Nuage de points
ggplot(data, aes(x=P, y=F))+geom_point(col="blue")+
geom_smooth(method = "lm")+
labs(title="Nombre de frere en fonction du poids", x="Poids",
y="Nombre de frere")
`geom_smooth()` using formula = 'y ~ x'
#Nuage de points
ggplot(data, aes(x=T, y=F))+geom_point(col="blue")+
geom_smooth(method = "lm")+
labs(title="Nombre de frere en fonction de la taille", x="Taille",
y="Nombre de frere")
`geom_smooth()` using formula = 'y ~ x'
#install.packages("GGally")
library(GGally)
Warning: le package 'GGally' a été compilé avec la version R 4.3.3
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
ggcorr(data[,c("P","T","F")])
table(data$S,data$C)
bleu brun gris noir vert
f 4 8 0 0 3
h 8 15 1 2 4
ggplot(data, aes(x=C))+geom_bar(aes(fill=S), position = "dodge")+
scale_fill_brewer(palette = "Pastel2")+
labs(title = "Distribution des individus suivant le sexe et la couleur des yeux",
x="Couleur des yeux", y="Frequence",fill="sexe")
library(grid)
library(vcd)
Warning: le package 'vcd' a été compilé avec la version R 4.3.3
#Diagramme en mosaic
mosaicplot(~ C + S ,
data = data,
color = TRUE,
las=1,
main = "Diagramme en mosaic du sexe et de la couleur des yeux")
#Matrice de contengence
ggplot(data, aes(y=S,x=C))+geom_bin_2d()+
labs(title = "Distribution des individus suivant le sexe et la couleur des yeux",
x="Couleur des yeux", y="sexe")
#Parametres sur les variables qualitatives
ass=assocstats(table(data$S, data$C))
print(ifelse(ass$cramer<=0.3,
paste("v cramer:",round(ass$cramer,2),"Le sexe et faiblement associé à la couleur des yeux", sep=" "),
ifelse(ass$cramer>0.3 & ass$cramer<=0.6,paste("v cramer:",round(ass$cramer,2),"Le sexe et moyennement associé à la couleur des yeux", sep = " "),
paste("v cramer:",round(ass$cramer,2),"Le sexe et moyennement associé à la couleur des yeux", sep = " "))))
[1] "v cramer: 0.2 Le sexe et faiblement associé à la couleur des yeux"
ggplot(data, aes(x=P, y=S))+geom_violin(aes(fill=S))+
labs(x = "Poids", y="Sexe", fill="Sexe",
title="Violinplot du poids des individus suivant le sexe")+
scale_fill_brewer(palette = "Pastel1")
ggplot(data, aes(y=P, x=C))+geom_boxplot(aes(fill=C),outlier.colour = "red",)+
labs(y = "Poids", x="Couleur des yeux", fill="Couleur des yeux",
title="Boxplot de la taille des individus suivant la couleur des yeux")+
scale_fill_brewer(palette = "Pastel2")
ggpairs(data[,c("P","T","F")])
ggally_density(data[,c("P","T","S")],aes(x=P, y=T, col=S))+labs(x="Poids", y="Taille", col="Sexe")
ggplot(data, aes(P, T))+stat_density_2d(aes(fill=S), geom = "polygon")+labs(x="Poids", y="Taille", col="Sexe")
On rejete H0 losque p-value < α (5% par defaut)
Pour tester l’homogeneite des variances (égalité des variances), de nombreux tests peuvent être utilisés notament :
Test F : Comparez les variances de deux groupes. Les données doivent être normalement distribuées.
Test de Bartlett : Comparer les variances de deux groupes ou plus. Les données doivent être normalement distribuées.
Le test de Levene : Une alternative robuste au test de Bartlett qui est moins sensible aux écarts de normalité.
Test de Fligner-Killeen : un test non paramétrique qui est très robuste contre les écarts de normalité.
NB : Il est à noter que le test de Levene est le plus couramment utilisé dans la littérature.
#install.packages("tidyverse")
#install.packages("conflicted")
#install.packages("rstatix")
library(conflicted)
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ lubridate 1.9.3 ✔ tibble 3.2.1
✔ purrr 1.0.2 ✔ tidyr 1.3.0
library(rstatix)
Warning: le package 'rstatix' a été compilé avec la version R 4.3.3
library(ggpubr)
#On visualise la normalité de la variable
ggqqplot(data, x="P")+labs(title="qqplot du poids de l'ensemble des individus")
ggqqplot(data, x="P", facet.by = "S")+labs(title="qqplot du poids de chaque sexe")
Les points étant presque tous situés autour de la droite de reference, nous pouvons supposer une normalité pour la variable poids. Cette normalité est vérifié par le test de shapiro wilk.
data %>% group_by(S) %>% shapiro_test(P)
NB : Si la taille de votre échantillon est supérieure à 50, le tracé QQ normal est préféré car avec des échantillons plus grands, le test de Shapiro-Wilk devient très sensible même à un écart mineur par rapport à la normalité.
var.test(P~S, data)
F test to compare two variances
data: P by S
F = 0.53037, num df = 14, denom df = 29, p-value = 0.2111
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.2251757 1.4533624
sample estimates:
ratio of variances
0.5303741
La P-Value > 5% indique que l’on accepte l’hypothèse H0 et donc que les variances des poids des deux sexes (h et f) sont égales.
Cette section décrit comment comparer plusieurs variances dans R à l’aide des tests de Bartlett, Levene ou Fligner-Killeen.
Hypothèses statistiques . Pour tous ces tests qui suivent, l’hypothèse nulle est que les variances de toutes les populations sont égales, l’hypothèse alternative est qu’au moins deux d’entre elles diffèrent. Par conséquent, des valeurs p inférieures à 0,05 suggèrent que les variances sont significativement différentes et que l’hypothèse d’homogénéité de la variance a été violée.
bartlett.test(P~S, data=data)
Bartlett test of homogeneity of variances
data: P by S
Bartlett's K-squared = 1.6963, df = 1, p-value = 0.1928
LeveneTest(P~S, data=data)
fligner.test(P~S, data=data)
Fligner-Killeen test of homogeneity of variances
data: P by S
Fligner-Killeen:med chi-squared = 0.20196, df = 1, p-value = 0.6531
Ce test permet de comparer les moyennes de deux groupes ou d’un groupe avec une valeur de reference connue
Dans cette situation nous avons un echantillon {x1, x2, …., xn} et on veut tester si la moyenne de cette echantillons est egale, superieur ou inferieur a une valeur connue.
On affirme que le poids moyen de la population est de 65kg. On utilise un t test pour verifier cette hypothese.
H0: mu=65
H1: mu<>65
mean(data$P)
[1] 65.08889
res<-t.test(data$P, mu = 65)
res
One Sample t-test
data: data$P
t = 0.066594, df = 44, p-value = 0.9472
alternative hypothesis: true mean is not equal to 65
95 percent confidence interval:
62.39882 67.77896
sample estimates:
mean of x
65.08889
On observe que la p-value > 5% donc on a bien le poids moyen de la population égale à 65 avec 5% de chance de se tromper.
res<-t_test(data, P~1, mu=65)
# Create the box-plot
bxp <- ggboxplot(
data$P, width = 0.5, add = c("mean", "jitter"),
ylab = "Weight (g)", xlab = FALSE
)
Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
Warning: The `fun.ymin` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun.min` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
Warning: The `fun.ymax` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun.max` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.
# Add significance levels
bxp + labs(subtitle = get_test_label(res, detailed = TRUE))
De même on peut tester si la taille moyenne de la population est superieure a 150cm
H0: mu=150
H1: mu > 150
mean(data$T)
[1] 174.3111
t.test(data$T, mu=150, alternative = "greater")
One Sample t-test
data: data$T
t = 18.456, df = 44, p-value < 2.2e-16
alternative hypothesis: true mean is greater than 150
95 percent confidence interval:
172.0978 Inf
sample estimates:
mean of x
174.3111
On observe que la p-value < 5%. Donc avec un risque de 5% de se tromper la taille de la population est superieur a 150cm
En utilisant cette fois le package rstatix on peut aussi tester pareillement si le nombre moyen de frere est inferieur a 3
H0: mu=3
H1: mu < 3
stat.test=data %>% t_test(F~1, mu=3, alternative = "less", detailed = TRUE)
stat.test
On observe ici qu’avec une confiance de 95% que le nombre de frère est bien inferieur à 3.
ggdensity(data, x = "F", rug = TRUE, fill = "lightgray") +
stat_central_tendency(type = "mean", color = "red", linetype = "dashed") +
geom_vline(xintercept = 3 , color = "blue", linetype = "dashed") +
labs(subtitle = get_test_label(stat.test, detailed = TRUE), x="Nombre de frere")
Le test t pour les écahntillons indépendants permet comme son nom l’indique de comparer les moyennes de deux échantillons biens différents. Il est aussi appelé t test pour un ecahntillon apparié.
Prérequis :
Une variable dependante catégorielle
Une variable indépendante numérique
Indépendance entre les observations de chaque catégorie
Données normalement distribuées
Vérifier l’égalité des variances (en fonction des résultats du test de levene mettre var.equal = TRUE/FALSE)
data %>% group_by(S) %>%
get_summary_stats(P, type = "mean_sd")
On peut supposer que le poids des hommes est superieur a celui des femme
H0: mu(h)=mu(f)
H1: mu > mu(f)
c’est a dire tester si
H0: mu(h)-mu(f)=0
H1: mu - mu(f)>0
#var.equal = TRUE car les variances sont égales d'après le test de levene
t.test(P~S, data=data, alternative="less",
var.equal = TRUE)
Two Sample t-test
data: P by S
t = -4.4418, df = 43, p-value = 3.074e-05
alternative hypothesis: true difference in means between group f and group h is less than 0
95 percent confidence interval:
-Inf -6.546832
sample estimates:
mean in group f mean in group h
58.06667 68.60000
Donc on peut affirmer avec 95% de chance qu’en moyenne les hommes pesent plus que les femmes
res<-t_test(data, P~S, alternative="less", var.equal = TRUE)
res <- res %>% add_xy_position(x = "S")
# Create the box-plot
bxp <- ggboxplot(
data, x="S", y="P", width = 0.5, add = c("mean"),
ylab = "Poids (kg)", xlab = "Sexe"
)
# Add significance levels
bxp +
stat_pvalue_manual(res, tip.length = 0) +
labs(subtitle = get_test_label(res, detailed = TRUE), title = "Boxplot du poids en fonction du sexe")
On peut évaluer cette différence en utilisant le d de cohens.
Ici on veut comparer les moyennes de deux mesures effectués sur des mêmes individus. Par exemple lorsqu’on veut comparer les mesures de tension arterielles des patients avant et après un traitement X.
#install.packages("datarium")
library(datarium)
Warning: le package 'datarium' a été compilé avec la version R 4.3.3
data("mice2", package = "datarium")
head(mice2, 3)
#Transformation des données en une colonne pour la classe et une pour les données mesurées
mice2.long <- mice2 %>%
gather(key = "group", value = "weight", before, after)
head(mice2.long)
#Resumé des données
mice2.long %>%
group_by(group) %>%
get_summary_stats(weight, type = "mean_sd")
#Execution du test
res<-t_test(weight~group, data=mice2.long, paired = TRUE, detailed = TRUE)
res
On peut conclure que le poids est significativement différent avant et après le traitement car la p-value est < 5%.
On peut évaluer cette différence en utilisant le d de cohens comme suit :
cohens_d(weight~group, data=mice2.long, paired = TRUE)
Donc l’écart entres les poids moyennes avant et après est d’environ 8.08
bxp <- ggpaired(mice2.long, x = "group", y = "weight",
order = c("before", "after"),
ylab = "Weight", xlab = "Groups")
res <- res %>% add_xy_position(x = "group")
bxp +
stat_pvalue_manual(res, tip.length = 0) +
labs(subtitle = get_test_label(res, detailed= TRUE))
Ce test d’applique a des variables binomiales telles que le sexe (homme et femme). Il permet de comparer une proportion d’un échantillon (exemple la proportion des hommes) avec une valeur de reference.
Exemple :
On suppose qu’il ya autant d’homme que de femme dans la population
H0 : P(h)=P(f)
H1 : P(h)<>P(f)
table(data$S)
f h
15 30
binom.test(table(data$S))
Exact binomial test
data: table(data$S)
number of successes = 15, number of trials = 45, p-value = 0.0357
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.2000057 0.4895036
sample estimates:
probability of success
0.3333333
p-value < 5% on rejete l’hypothèse nulle et donc la proportion des hommes est significatvement différente de celle des femmes
Le test du Khi-deux est un test d’hypothèse utilisé pour déterminer s’il existe une relation entre deux variables catégorielles.
data$PR<-cut(data$P, breaks=5)
ggplot(data, aes(x=PR, fill=S))+geom_bar(position = "dodge")
On pourrait supposer une liaison entre le sexe et la tranche de poids
khi_test<-chisq.test(x=data$S, y=data$PR)
Warning in chisq.test(x = data$S, y = data$PR): L’approximation du Chi-2 est
peut-être incorrecte
print("Distribution observe")
[1] "Distribution observe"
khi_test$observed
data$PR
data$S (46.9,57.2] (57.2,67.4] (67.4,77.6] (77.6,87.8] (87.8,98.1]
f 6 8 1 0 0
h 2 11 15 1 1
print("Distribution attendue")
[1] "Distribution attendue"
khi_test$expected
data$PR
data$S (46.9,57.2] (57.2,67.4] (67.4,77.6] (77.6,87.8] (87.8,98.1]
f 2.666667 6.333333 5.333333 0.3333333 0.3333333
h 5.333333 12.666667 10.666667 0.6666667 0.6666667
khi_test
Pearson's Chi-squared test
data: data$S and data$PR
X-squared = 13.189, df = 4, p-value = 0.01039
NB : Les distributions théoriques (attendue) doivent toutes être supérieures à 5
Le test U de Mann-Whitney peut être utilisé pour vérifier s’il existe une différence entre deux échantillons (groupes), et les données ne doivent pas nécessairement être distribuées normalement. Pour pouvoir calculer un test U de Mann-Whitney, il faut disposer de deux échantillons aléatoires indépendants présentant au moins des caractéristiques à échelle ordinale.
H0 : Il n’existe pas de différence entre les groupes (du point de vu des tendances centrales)
group_by(data,S) %>%
summarise(
count = n(),
median = median(P, na.rm = TRUE),
IQR = IQR(P, na.rm = TRUE))
wilcox.test(P~S, data = data)
Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): impossible
de calculer la p-value exacte avec des ex-aequos
Wilcoxon rank sum test with continuity correction
data: P by S
W = 56, p-value = 4.739e-05
alternative hypothesis: true location shift is not equal to 0
Le test de Wilcoxon (test du rang de signe de Wilcoxon) vérifie si les valeurs moyennes de deux groupes dépendants diffèrent significativement l’une de l’autre.
Le test de Wilcoxon est un test non paramétrique et est donc soumis à beaucoup moins de conditions prélables que son homologue paramétrique, le test t pour échantillons dépendants. Par conséquent, dès que les conditions limites du test t pour échantillons dépendants ne sont plus remplies, le test de Wilcoxon est utilisé.
Conditions du test de Wilcoxon
Le test de Wilcoxon étant un test non paramétrique, il n’est pas nécessaire que les données soient normalement distribuées. Toutefois, pour calculer un test de Wilcoxon, les échantillons doivent être dépendants. Les échantillons dépendants sont présents, par exemple, lorsque les données sont obtenues à partir de mesures répétées ou lorsqu’il s’agit de paires dites naturelles.
Mesure répétée : une caractéristique d’une personne, par exemple son poids, a été mesurée à deux moments différents.
Couples naturels : les valeurs ne doivent pas nécessairement provenir de la même personne, mais de personnes qui vont ensemble, par exemple avocat/client, épouse/mari et psychologue/patient. Bien entendu, il ne s’agit pas nécessairement de personnes.
Indépendance : le test de Wilcoxon suppose l’indépendance, c’est-à-dire que les observations appariées sont tirées au hasard et de manière indépendante.
En outre, la forme de la distribution des différences entre les deux échantillons dépendants doit être approximativement symétrique.
Si les données ne sont pas disponibles par paires, le test U de Mann-Whitney est utilisé à la place du test de Wilcoxon.
# Load the data from datarium package
data("genderweight", package = "datarium")
# Show a sample of the data by group
set.seed(123)
genderweight %>% sample_n_by(group, size = 2)
genderweight %>%
group_by(group) %>%
get_summary_stats(weight, type = "median_iqr")
bxp <- ggboxplot(
genderweight, x = "group", y = "weight",
ylab = "Weight", xlab = "Groups", add = c("jitter")
)
bxp
stat.test <- genderweight %>%
wilcox_test(weight ~ group) %>%
add_significance()
stat.test
stat.test <- stat.test %>% add_xy_position(x = "group")
bxp +
stat_pvalue_manual(stat.test, tip.length = 0) +
labs(subtitle = get_test_label(stat.test, detailed = TRUE))
Le test ANOVA (ou Analyse de Variance en francais) est utilisé pour comparer la moyenne de plusieurs groupes. Le terme ANOVA est un peu trompeur. Bien que le nom de la technique fasse référence aux variances, l’objectif principal de l’ANOVA est d’étudier les différences de moyennes.
Une extension du test t pour échantillons indépendants pour comparer les moyennes dans une situation où il y a plus de deux groupes. Il s’agit du cas le plus simple de test ANOVA où les données sont organisées en plusieurs groupes selon une seule variable de regroupement (également appelée variable factorielle). D’autres synonymes sont : ANOVA à 1 voie , ANOVA à un facteur et ANOVA inter-sujets
data("PlantGrowth")
PlantGrowth <- PlantGrowth %>%
reorder_levels(group, order = c("ctrl", "trt1", "trt2"))
PlantGrowth %>%
group_by(group) %>%
get_summary_stats(weight, type = "mean_sd")
ggboxplot(PlantGrowth, x = "group", y = "weight")
Les valeurs aberrantes peuvent être facilement identifiées à l’aide
des méthodes de diagramme en boîte, implémentées dans la fonction
R identify_outliers()du package
rstatix.
PlantGrowth %>%
group_by(group) %>%
identify_outliers(weight)
Vérification de l’hypothèse de normalité.
Elle peut être fait soit :
A partir des tests de normalité dans chaque groupe
Ou a partir des résidus du modèle ANOVA
# Construction du modèle linéaire
model <- lm(weight ~ group, data = PlantGrowth)
# QQ Plot
ggqqplot(residuals(model))
Les valeurs sont proches de la droite de reférence, cela suppose une normalité des résidus
H0 : La variable est normalement distribué
#Test de normalité de shapiro wilk
shapiro_test(residuals(model))
Dans le cas de la normalité de chaque groupe on aurait
ggqqplot(data=PlantGrowth, x="weight", facet.by = "group")
PlantGrowth %>% group_by(group) %>%
shapiro_test(weight)
La P-value ici est > 5% donc le poids suit une distribution normale.
res.aov <-PlantGrowth %>% anova_test(weight ~ group)
Une ANOVA unidirectionnelle significative est généralement suivie de
tests post-hoc de Tukey pour effectuer plusieurs comparaisons par paires
entre les groupes. Fonction de la touche R :
tukey_hsd()[rstatix].
pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)
#Visualisation
pwc <- pwc %>% add_xy_position(x = "group")
ggboxplot(PlantGrowth, x = "group", y = "weight") +
stat_pvalue_manual(pwc, hide.ns = TRUE) +
labs(
subtitle = get_test_label(res.aov, detailed = TRUE),
caption = get_pwc_label(pwc)
)
NB : Le test ANOVA unidirectionnel classique nécessite une hypothèse de variances égales pour tous les groupes. Dans le cas contraire, il est conseillé d’utiliser :
Le test unidirectionnel de Welch est une alternative à l’ANOVA unidirectionnelle standard dans les situations où l’homogénéité de la variance ne peut pas être supposée (c’est-à-dire que le test de Levene est significatif).
Dans ce cas, le test post hoc de Games-Howell ou les tests t par paires (sans hypothèse de variances égales) peuvent être utilisés pour comparer toutes les combinaisons possibles de différences de groupe.
#Dans ce cas voici les fonctions à utiliser
#welch_anova_test()
#games_howell_test()
Utilisée pour évaluer simultanément l’effet de deux variables de regroupement différentes sur une variable de résultat continue. D’autres synonymes sont : deux plans factoriels , anova factorielle ou ANOVA bidirectionnelle entre sujets .
data("jobsatisfaction", package = "datarium")
jobsatisfaction %>%
group_by(gender, education_level) %>%
get_summary_stats(score, type = "mean_sd")
ggplot(jobsatisfaction)+geom_boxplot(aes(x=gender, y=score, fill=education_level))
model <- lm(score ~ gender*education_level,
data = jobsatisfaction)
ggqqplot(residuals(model))
shapiro_test(residuals(model))
jobsatisfaction %>%
group_by(gender, education_level) %>%
shapiro_test(score)
ggqqplot(jobsatisfaction, "score", ggtheme = theme_bw()) +
facet_grid(gender ~ education_level)
jobsatisfaction %>% levene_test(score ~ gender*education_level)
res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov
Une interaction bidirectionnelle significative indique que l’impact d’un facteur (par exemple, le niveau_d’éducation) sur la variable de résultat (par exemple, le score de satisfaction au travail) dépend du niveau de l’autre facteur (par exemple, le sexe) (et vice versa). Ainsi, vous pouvez décomposer une interaction bidirectionnelle significative en :
Effet principal simple : exécuter un modèle unidirectionnel de la première variable à chaque niveau de la deuxième variable,
Comparaisons simples par paires : si l’effet principal simple est significatif, effectuez plusieurs comparaisons par paires pour déterminer quels groupes sont différents.
Pour une interaction bidirectionnelle non significative , vous devez déterminer si vous avez des effets principaux statistiquement significatifs à partir du résultat de l’ANOVA. Un effet principal significatif peut être suivi par des comparaisons par paires entre les groupes.
model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
group_by(gender) %>%
anova_test(score ~ education_level, error = model)
L’effet principal simple du « niveau_d’éducation » sur le score de satisfaction au travail était statistiquement significatif tant pour les hommes que pour les femmes (p < 0,0001).
En d’autres termes, il existe une différence statistiquement significative dans le score moyen de satisfaction au travail entre les hommes ayant fait des études scolaires, collégiales ou universitaires, F(2, 52) = 132, p < 0,0001. La même conclusion est vraie pour les femmes , F(2, 52) = 62,8, p < 0,0001.
NB : Il convient de noter que la signification statistique des analyses simples de l’effet principal a été acceptée à un niveau alpha ajusté par Bonferroni de 0,025. Cela correspond au niveau actuel auquel vous déclarez la signification statistique (c’est-à-dire p < 0,05) divisé par le nombre d’effets principaux simples que vous calculez (c’est-à-dire 2).
https://datatab.fr/tutorial/get-started
https://3mmarand.github.io/comp4biosci/
https://epirhandbook.com/en/index.html
https://openintro-ims.netlify.app/
https://larmarange.github.io/analyse-R/graphiques-bivaries-ggplot2.html
https://www.datanovia.com/en/lessons/how-to-do-a-t-test-in-r-calculation-and-reporting/
https://www.umr-lastig.fr/paul-chapron/resources/cours_site/galerie-de-graphiques-avec-ggplot.html
https://www.datanovia.com/en/fr/lessons/test-dhomogeneite-des-variances-dans-r/
https://jokergoo.github.io/ComplexHeatmap-reference/book/introduction.html
http://www.sthda.com/english/wiki/descriptive-statistics-and-graphics
http://www.sthda.com/english/wiki/ggpubr-create-easily-publication-ready-plots